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We provide a comprehensive presentation of the Hierarchical Reference Theory (HRT) in the 
smooth cut-off formulation. A simple and self-consistent derivation of the hierarchy of differential 
equations is supplemented by a comparison with the known sharp cut-off HRT. Then, the theory is 
applied to a hard core Yukawa fluid (HCYF): a closure, based on a mean spherical approximation 
ansatz, is studied in detail and its intriguing relationship to the self consistent Ornstein-Zernike 
" approximation is discussed. The asymptotic properties, close to the critical point are investigated 

. and compared to the renormalization group results both above and below the critical temperature. 

■ The HRT free energy is always a convex function of the density, leading to flat isotherms in the two- 

| phase region with a finite compressibility at coexistence. This makes HRT the sole liquid-state theory 

able to obtain directly fluid-fluid phase equilibrium without resorting to the Maxwell construction. 
The way the mean field free energy is modified due to the inclusion of density fluctuations suggests 
how to identify the spinodal curve. Thermodynamic properties and correlation functions of the 
HCYF are investigated for three values of the inverse Yukawa range: z = 1.8, z = 4 and z = 7 
where Monte Carlo simulations are available. The stability of the liquid-vapor critical point with 
respect to freezing is also studied. 
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PACS numbers: 64.60.F-, 61.20.Gy, 64.60.A-, 05.70.Fh 



I. INTRODUCTION 



Our understanding of the liquid- vapor phase transition was boosted, in the seventies, by the introduction in a statis- 
[ tical physics framework, of genuinely field theoretical methods: long wavelength effective actions and Renormalization 
£H . Group (RG) techniques The microscopic justification of these coarse grained approaches has been the subject of 
several studies in the eighties 0, 0, 0] paving the way to quantitative investigations of the phase diagrams of simple 
O . liquids and mixtures [f|. 

In particular, the Hierachical Reference Theory of fluids (HRT) [4J, |fj] can be considered as a successful attempt at 
developing the RG methods in the framework of a liquid-state theory. HRT is based on the microscopic hamiltonian of 
the fluid, but its RG-like structure endowes it with a number of features that are not shared by other theories of liquids: 
specifically, a non-trivial description of critical behavior that incorporates universality and scaling and, even more 
remarkably, a treatment of first-order phase transitions that does not violate the convexity of the free energy required 
by thermodynamic stability, thereby yielding rigorously flat pressure vs. density isotherms whenever different phases 
coexist at equilibrium. The region of phase coexistence is therefore straightforwardly obtained by HRT as the locus 
of diverging compressibility, and does not have to be recovered a posteriori by enforcing thermodynamic equilibrium 
via the Maxwell construction, that is instead necessary in other liquid-state theories, and may turn out to be quite 
hard to implement, especially in the critical region or in the case of binary fluids. 

Notwithstanding these qualities, the original formulation of HRT leaves room for improvement. Specifically, its most 
serious deficiency consists in the fact that it fails to predict the discontinuity of the inverse compressibility expected 
in systems with a scalar order parameter, when the coexistence curve is crossed at any noncritical temperature. In 
fact, HRT predicts a diverging compressibility not only inside the coexistence region, but also on its boundary, so 
that the coexistence and spinodal curves coincide 0]. Besides preventing one from recovering the correct behavior 
of the compressibility at coexistence, this deficiency also forbids the possibility for the system to survive in a single 
phase inside the coexistence region as a metastable state — an occurrence that would be allowed to states located 
between the coexistence and the spinodal. One might think of this as a minor problem, since the very concept of 
a spinodal curve has little to do with equilibrium thermodynamics. However, the instance of metastability and its 
breakdown via the process of spinodal decomposition are experimentally relevant, especially in the context of complex 
fluids. These systems, often characterized by short range, tunable potentials and/or competing interactions, now offer 
the possibility to tailor phase diagrams and thermophysical properties 0] giving rise to phase boundaries of different 
topologies. Non-equilibrium states play a major role in these systems, being rather ubiquitous in colloidal physics. 
In particular, the relevance of spinodal decomposition for the occurrence of attractive gel phases has recently been 
pointed out 0]. 

Here we present a new formulation of HRT, the smooth cut-off HRT, that is free from the aforementioned short- 
coming. We will detail the physical motivations at the basis of this generalization, together with the technical 
implementation of the smooth cut-off HRT to a rather flexible class of models: the Yukawa fluids. For this class of 
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potentials, the smooth cut-off HRT allows also to implement exactly the requirement that the particles be mutually 
impenetrable because of the hard-core part of the interaction, usually referred to as the core condition. A numerical 
study of the resulting equations allows to test the accuracy of the theory via the comparison with available numerical 
simulations. 

The concept at the heart of HRT is the close relationship between the amplitude of density fluctuations of given 
wavevector k and the corresponding Fourier component of the attractive part of the potential. This observation 
has a long history, dating back to the early works on the Random Phase Approximation (To| and to the Hubbard- 
Schofield mapping of a fluid model onto a scalar field theory [3j. In order to make this correspondence explicit, it is 
customary to split the two body inter-particle potential into the sum of a reference (mainly repulsive) part v R (r) and 
a residual attractive tail w(r). The formal perturbation series in powers of w(k) (hereafter a "tilde" denotes Fourier 
transforms) can be reinterpreted, at all orders, as the diagrammatic expansion of a suitable effective scalar field 
theory governing density fluctuations, thereby providing a bridge between standard liquid state theory and the field 
theoretical implementation of RG @ . According to the momentum space RG philosophy , the order parameter 
fluctuations on different length-scales arc gradually introduced into the system in an iterative way, starting from short 
wavelengths (i.e. large wavevectors) down to k = 0. As a consequence, the liquid-vapor phase transition, whose order 
parameter is just the k = component of density fluctuations pk, is inhibited at any intermediate step of the RG flow 
and is recovered only at the end of the iterative procedure. This sort of "regularization" implicit in the RG allows the 
use of rather crude approximations at each step of integration, still preserving the scaling properties of the theory. 
This procedure can be conveniently exported into liquid state theory by taking advantage of the previously outlined 
correspondence between the amplitude of density fluctuations pk and Fourier components of the potential w(k): the 
effects of w(k) have to be introduced into the model selectively in the wavevector k. Such a criterion, however, leaves 
an ample range of possibilities in the practical implementation of the RG procedure. The most natural choice is 
to define a sequence of fictitious interactions WQ(k) with vanishing Fourier components for k < Q. By following 
the flow of the physical properties of the model fluid as Q varies from infinity to zero, we mimic the action of the 
momentum space RG approach: this program has been pursued in the usual "sharp cut-off" formulation of HRT 
[3, @. However, in this way the two physically sensible models, the reference system (recovered for Q = oo) and the 
fully interacting fluid (corresponding to Q = 0) are connected by a sequence of artificial models, whose interaction 
is characterized by long range oscillatory tails due to the sharp discontinuity of WQ(k) in momentum space. Such 
a pathological interaction may pose some specific problem for the formulation of the approximations necessary to 
evaluate the physical properties of the model along the flow, i.e. at intermediate Q's. It would be desirable to define a 
different sequence of intermediate systems which, still maintaining the selective inclusion of wavevectors, correspond 
to finite range, regular interactions at each step. A possible solution of the problem, the smooth cut-off HRT, was 
already put forward several years ago [l2| but then the analysis did not go beyond a preliminary investigation on the 
critical properties predicted by this approach, later supplemented by a study of the first order transition [l3|. In the 
following we will review the formulation of the smooth cut-off HRT equations providing full details on the theory and 
a thorough discussion of the results, some of which have been already anticipated in Ref. [14| . 



II. THE SMOOTH CUT-OFF HRT EQUATIONS 



Following the program previously outlined, we consider a gas of classical particles in dimension d, interacting via a 
spherically symmetric two body potential v(r) written as the sum of a reference part and a tail w(r): 

v(r) = v R (r) + w(r) (I) 

Then we define a sequence of intermediate systems interacting with potential Vt(r) parametrized by t £ (0, oo) which 
interpolates between va(r) (for t = 0) and v(r) (as t — > oo). This procedure is able to capture the RG spirit whenever 
v t (r) effectively suppresses the Fourier components of w(r) with k < Q 4- e _t . This effect can be obtained by different 
definitions of vt (r) . Here we concentrate on the following class of parametrizations fl2| : 

Vt(r) = v R (r) + w t (r) 

= v R (r) + [w(r) - t()(t) e- dt w(r e-*)] (2) 

where tf)(t) is a monotonically decreasing switching function, to be specified later, vanishing as t — > oo and satisfying 
■0(0) = 1. The Fourier transform of Wt(r) is given by: 

w t {k)=w{k)-i){t)w{ke t ) (3) 

For a short range interaction characterized by an attractive tail with Fourier components w(k), negligible beyond 
a certain wave vector Q*, the intermediate potential w t (k) coincides with w(k) for k > Q* e~* while ip(t) < I 
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inhibits the Fourier components at k < Q* . Therefore the definition ((2|) does satisfy the requirements discussed in 
the introduction, gives rise to a smooth, regular interaction at all t's and may represent the starting point for a 
microscopic implementation of the RG approach. 

First order perturbation theory (l5j provides the exact expression for the change in the free energy density of this 
model when t is varied: 

dA t a p 2 f , , . dw t {r) 

where A t is minus the excess free energy density divided by ksT, gt{r) is the radial distribution function of the 
intermediate system and (3 = 1/fcgT. This equation relates the change in the free energy with the change in the 
interaction but requires the knowledge of the two body correlations. An analogous equation for the pair correlations 
can be obtained starting from the formal diagrammatic expansion of the free energy functional in powers of the 
interaction, as derived in Ref. 



ic| 4) (k, -k, q, -q) + cf (k, -q, q k) (|k - q|) cf } (-k, q, k q) 



dt 
(5) 

where Cj (ki . . . k„) is the n-particle direct correlation function (including the ideal gas contribution)defined in terms 
of the n th functional density derivative of the total free energy of the partially interacting system |6j. In particular, 

(2) 

c t (k) is related to the usual direct correlation function ct(fc) by 

~cf\k) = --+c t {k). 
P 

(2.) ~ (2.) 

The Ornstein-Zernike equation connects c\ (k) to F£ (fc), which equals the structure factor S t (k) multiplied by p: 

Fl 2 \k) = pS t (k) = - P -— = p + p 2 f dre^[g t (r)-l}. 

l-pc t (fc) J 

Equations (|4l5p are the first two terms of an exact hierarchy of differential equations for the free energy and the many 
body direct correlation functions of the model. Together with the definition of Wt(r) $2$ this hierarchy defines the 
smooth cut-off formulation of HRT. It is interesting to notice that when Eq. ([5]) is evaluated at vanishing external 
momentum (fc = 0), it becomes equivalent to the second density derivative of the HRT equation for the free energy 
(j4|). An analogous correspondence takes place throughout the full hierarchy, the (n + l) th differential equation being 
related to the density derivative of the n . This coincidence is not accidental, because the direct correlation functions 
are just functional derivatives of the free energy, and the functional derivative evaluated at zero momentum reduces 
to the partial derivative with respect to the uniform density. This shows that the first HRT equation (j4j does contain 
information on the structure of the full hierarchy by virtue of the exact relationship between the free energy and the 
direct correlation function implied by the compressibility sum rule: 

In order to exploit the connection of the HRT equations with the RG theory, it is convenient to rescale appropriately 
the set of direct correlation functions. As previously discussed, the switching parameter t provides the wavevector 
scale of density fluctuations Q-r- e~*. At criticality, according to general scaling arguments the expected behavior 
of the n-particle direct correlation function (also known as one particle irreducible n-point function) is then given by 

C 1 n) ( qi . . . q„) = _ e [^- 2+ "M 4 ^(qie* • ■ ■ q^') (7) 

where rj is the critical exponent related to the anomalous dimension The scaling form which characterizes 

correlations at the critical point is recovered provided the function u[ n \x.i . ..x„) tends to a finite limit as t — > oo. 
By substituting Eq. ([7]) into the full (infinite) set of HRT equations, it is easy to show that the explicit dependence 

(n) 

on the parameter t disappears from the equations for the scaling functions u t if the switching function is given by 
4>(t) — e~( 2 ~ ?? )*. With this choice of ip(t), the full hierarchy, rescaled through the definition {7J, admits a fixed point 
solution providing the formal link with RG methods and allowing for a microscopic derivation of the known results on 
scaling laws and critical exponents fl2| in full agreement with the known dimensionality expansion [l[ . Note however 
that the scaling form is expected to hold only at small wavevectors, i.e. at large t, implying that the HRT equations 
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provide a microscopic implementation of the RG procedure for any choice of ip(t) characterized by the asymptotic 
behavior: 

Km e (2 -*V(*) = ^oo (8) 

t — >oo 

with finite and non zero. 

A. Sharp vs. smooth cut-off procedure 

The original implementation of HRT |4jwas based on a different switching-on of the attractive interaction, inspired 
by the momentum shell integration RG (llj . The sequence of intermediate systems are defined by an attractive part 
of the potential wq(t) whose Fourier components are set to zero below a cut-off wave vector Q. By varying Q from 
infinity to zero vq(t) = Vr(t) + wq(t) interpolates between the reference system and the fully interacting one. In this 
case, however, the "slice" of potential Swq = wq_<jq — wq included at each step is characterized by a finite Fourier 
transform w(Q) in an infinitesimal domain (Q — dQ, Q). Such a singular form of Swq has profound consequences on 
the structure of the HRT equations: now, first order perturbation theory is not sufficient to evaluate the change in 
the free energy when Q is infinitesimally decreased and a resummation of the ring diagrams is required [4j. Therefore, 
the HRT equations in the sharp cut-off formulation do not coincide with Eqs. (|4I5[) : in particular, the momentum 
integrations are limited to the shell q — Q suggesting a decoupling of the long wavelength properties of the fluid model. 
In the critical region, the sharp cut-off HRT hierarchy acquires a form independent of the interparticle interaction, 
leading in a natural way to the concept of universality of the critical behavior. This appealing feature is not share d by 
the smooth cut-off hierarchy where the proof of universality requires a more elaborate analysis of the equations [12| . 
On the other hand, the intermediate potentials wq (r) in the sharp cut-off formulation display unphysical long range 
oscillatory tails as a consequence of the discontinuity in Fourier space. This circumstance brings about some difficulty 
in the formulation of an approximate closure to the infinite hierarchy of differential equations, due to the absence of 
an underlying physical intuition on the properties of fluid models characterized by such a class of potentials. The 
simple Ornstein-Zernike (OZ) closure to the first equation of the hierarchy adopted in all the previous investigations 
of HRT, although provides non classical critical exponents and scaling law together with flat isotherms at coexistence, 
has been proved to induce some artificial feature in the description of first order phase transitions leading to the 
divergence of the compressibility at the phase boundary, i.e to the coincidence of the spinodal with the binodal line 
Here we will show how the same class of OZ closures instead provides a satisfactory description of the phase 
diagram of simple fluids within the smooth cut-off formulation of the HRT equation for the free energy (|2l4l8p . 



III. A SIMPLE CLOSURE FOR YUKAWA POTENTIALS 



The OZ closure we are going to study (named HRT-OZ) amounts to parametrize the radial distribution function 
for the intermediate system gt(r) appearing in the first equation of the HRT hierarchy (fj[|). Here we will focus on a 
simple yet flexible class of interactions: the hard core plus Yukawa potentials defined as the sum of a pure hard core 
term of diameter a and an attractive Yukawa tail of inverse range z 

w(r) = -e (9) 

r 

In the following, a and e/fc^ will be taken as units of length and temperature respectively (i.e. we set a = 1, e = fcg). 

As previously noticed, in order to ensure that a memory of the full structure of the HRT hierarchy is contained 
in the free energy equation f4|, the compressibility sum rule ([6]) has to be satisfied within the parametrization. A 
simple way to implement such a requirement is to adopt an approximate form of the correlation functions coming 
from some accurate liquid state theory suitably generalized in order to allow for the consistency between free energy 
and correlations ([6]). In order to take advantage of the analytical form available for the correlation functions of a 
Yukawa fluid, we have chosen to use, as a supporting liquid state theory, the Mean Spherical Approximation (MSA) 
defined by the usual core condition (.gt(r) = for r < 1) supplemented by the explicit form of the direct correlation 
function outside the core: 

c*(r) = —/3(w t (r) +\ t w(r)) 

= -f3(l + \ t )w(r) +/3V>(t)e~ dt w(re"*) for r>l (10) 

where A t is implicitly defined by the compressibility sum rule (d|) in terms of the excess free energy A t . The radial 
distribution function and the direct correlation function are connected via the OZ equation. Within MSA the direct 
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correlation function is always analytic in k 2 at small k, even at the critical point, implying that the critical exponent 
77 vanishes within this class of approximations. For the HCYF, Eq. (fTU|) implies a direct correlation function written 
as the sum of two Yukawa's of different ranges z% and amplitudes Ki'. 

p— 21(1 — 1) e~ Z2 (' — ^ 
c t {r)=K l + K 2 (11) 



outside the particle core (for i.e. r > 1), with 



Z\ = z 

z 2 = ze~* (12) 



K x = 



1 + A 



T 

e ^(i-e-*) 



X 2 = -—^e-^ (13) 

In three dimensions, the MSA integral equation can be solved analytically for this class of parametrizations and the 
solution reduces to a set of algebraic equations |l6[. Within HRT-OZ, this set is coupled to the evolution equation 
(|4"|) via the compressibility sum rule ^ which implicitly defines the parameter X t present in K\. The explicit form 
of the equations, together with a detailed description of the numerical scheme adopted for their resolution, can be 
found in Appendix A. Here we remark that, as discussed in detail in Appendix A, the parametrization (|10[) also 
allows to improve the performance of the MSA (and then of the HRT-OZ) in the high density region, where the hard 
sphere direct correlation function is known to display a non negligible tail outside the particle core. If this additional 
contribution to Ct(r) is represented as a Yukawa tail of the same inverse range z as that of the attractive part of the 
potential, its presence is equivalent to a redefinition of the parameter A* for t = 0, whose value is determined by Eq. 
([6]) if the Carnahan-Starling expression for the hard sphere free energy is used. 

A. HRT-OZ vs. SCOZA 

The Self consistent Ornstein Zernike Approximation (SCOZA) is a well known liquid state theory which proved 
successful in the determination of phase diagrams in Yukawa fluids 17]. Originally, SCOZA has been formulated as 
a generalization of the MSA expression for the direct correlation function outside the core: 

c(r) = -/3\w{r) (14) 

where the parameter A is determined by requiring thermodynamic consistency between the internal energy and the 
compressibility routes to thermodynamics. In terms of the excess free energy density (divided by — fc^T) A, this 
condition reads: 

-QP = - Y J dr 9{r)w{r) 
B 2 A f 

= (15) 

The two exact expressions (fl"5|) . together with the core condition, g{r) = for r < 1, and the parametrization (fTi]) 
for r > 1, give rise to a partial differential equation for the free energy as a function of p and (3. The close similarity 
between Eqs. (fT5|) and (|4|6p suggests that SCOZA can be somehow related to the smooth cut-off formulation of HRT. 
In fact, by defining the special sequence of intermediate potentials Wt(r) as: 

Wt (r) = t w(r) (16) 

the HRT-OZ equations (|4I6|) exactly reproduce SCOZA, showing that the two theories differ just by the switching- 
on procedure of the attractive pat of the interaction: while SCOZA tunes the amplitude of w(r), HRT changes 
simultaneously both amplitude and range of the potential. The long range repulsive tail introduced by HRT during 
the integration, allows to inhibit phase separation at all finite £'s thereby mimicking the RG approach. 

In order to emphasize the key role played by the weak long range repulsive tail added to the physical interaction in 
the smooth cut-off HRT, let us consider the simple MSA for a Yukawa fluid, which, at low temperatures, is known to 
display a region inside the spinodal where no solution exists. However, as shown in Fig. [TJ if the physical attractive 
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FIG. 1: Dimensionless inverse compressibility St(0) 1 of a z — 1.8 Yukawa fluid within MSA below the critical temperature 
(T = 1). The potential is given by © with ip(t) = (cosht)" 2 (see Eq. |63j at t = 2 (dotted line) and t = 10 (full line). Vertical 
dashed lines bracket the coexistence region: here the t = 10 isotherm is smaller than 10~ r . The standard MSA result for a 
single Yukawa coincides with the full line outside coexistence while no solution exists inside the coexistence region. 

Yukawa potential is supplemented by a long range repulsive Yukawa contribution of the form |2|). the spinodal curve 
disappears and a solution of the MSA equation exists in the whole phase diagram. When t — > oo, the repulsive term 
in the potential becomes vanishingly small: the isotherm becomes flat and the compressibility diverges in the whole 
region inside the spinodal of the single Yukawa, thereby enforcing the convexity of the free energy. Therefore, cutting 
off long wave length density fluctuations by the addition of a small, long ranged repulsive tail appears to be an easy 
way to implement one of the most important features of equilibrium statistical mechanics within an approximate 
liquid state theory. Note however that within MSA, this kind of "regularization" of the physical two body interaction 
leads to the collapse of the spinodal and coexistence curve. 

On the basis of these observations, the smooth cut-off HRT-OZ approach is then expected to provide a convex 
free energy: in the following Section we will show that HRT-OZ enjoys other relevant physical properties of first and 
second order phase transitions. 



As already noticed, the critical properties of the theory are determined by long wavelength fluctuations which come 
into play at the last stages of the i-evolution, i.e. in the limit t — > oo. The full HRT-OZ equations, in the critical 
region and at large i's considerably simplify: remarkably, the formal structure of the asymptotic HRT-OZ equation 
is in fact identical to that obtained by a much simpler closure to the first equation of the hierarchy (HJ) based on the 
Random Phase Approximation (RPA), i.e. neglecting the core condition. In fact, according to the critical phenomena 
paradigm, short range correlations should not affect the physics at long wavelengths which dominates in the critical 
region. And this property is explicitly satisfied within HRT-OZ. 

In the simplified RPA form, the direct correlation function reads 



for all r's. As usual, the parameter X t is determined by the compressibility sum rule ©■ Appendix B contains the 
details of the asymptotic analysis of (jl]) under both closures (fT7| and (fT0|) . Here we simply discuss the physical 
contents of the result. It is convenient to first define a modified free energy density (in units of fc^T) At which 
encompasses both the ideal gas term (A ld ) and the mean field contribution due to that part of the interaction not 
already included in A t : 



IV. ASYMPTOTIC ANALYSIS 



ct{r) = c R (r) - (3 (w t (r) + A t w(r)) 



(17) 



At = A t + A id 



(18) 



The evolution equation ([¥]) under the closure Q17p for a three dimensional Yukawa fluid in the t — > oo limit and in the 



7 



critical region becomes (see Appendix B): 

_3 

-1/2. 



dA t = _ 3 t£l 

at e 4tt 



p 



(19) 



where 



* 92 At (20) 



is a renormalized inverse compressibility and the non universal parameter p < 1 is the ratio between the curvature 
of the potential and of the direct correlation function divided by ipoo ©, as shown in Eqs. (|66l78p . It is convenient 
to remove an analytic contribution and a regular pre-factor in Eq. (|19p by a suitable rescaling of the free energy 
At — > f&t and of the density p — > ip as detailed in Eq. ([68)1 : 

/ z 5 \ 1/2 

p - pc = ip \WfNP) ■ (21) 

The resulting evolution equation reads: 



where (f^U|) is now expressed as 



.2* ^ 



Notice that, according to the definition (|2H)l . the variable x may attain a finite limit (x*) for t — ► oo only if the 
compressibility diverges, meaning that the fully interacting system is either at the critical point or in the two phase 
region. The asymptotic Eq. (|22[) has been numerically integrated for an initial condition at t = of the Landau- 
Ginzburg type, i.e. appropriate for a self-interacting scalar field theory: 

^o((f) — *o(0) + rip 2 +up 4 (24) 

Besides p, which retains information on the range of the potential, this effective model is defined by two parameters: r 
is a measure of temperature (shifted by the mean field critical value) while u > is assumed to be state independent. 
The typical behavior of the rescaled inverse compressibility x(<p, r) during the "time" evolution determined by Eq. 
(|22|) is shown in Fig. [2] for p = 0.5, u = 0.1 and three states on the critical isochore ip = 0: above, below and very 
close to the critical temperature r c . Above r c the divergence of x shows that the inverse compressibility attains a 
finite, positive value, as expected in supercritical fluids. Close but still above the critical point, the evolution of x 
displays a plateau at x* < before diverging again. Instead, for r < r c the evolution crosses x* and asymptotically 
sets to a finite, negative value x< implying that the inverse compressibility vanishes. Inside the coexistence curve the 
fluid is therefore characterized by flat isotherms, i.e. the Maxwell construction is recovered via the inclusion of long 
wavelength fluctuations. Such a behavior of the "time" evolution suggests the possible presence of two distinct fixed 
points of the evolution equation, x* and x<, governing the critical and subcritical region respectively (while the "high 
temperature fixed point" is located at x = +oo). 

A. Critical behavior 

The singular behavior of the thermodynamic properties in the critical region can be obtained from the asymptotic 
form of the evolution equation (|22p following the fixed point analysis developed in the RG theory. As a first step 
we eliminate from Eq. (|22[) the explicit dependence on the evolution parameter t by a i-dependent rescaling of the 
density and free energy density: 

C = e?<p 

r t (c) = ~ *t(Q)] e 3t (25) 
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FIG. 2: Evolution of the scaled inverse compressibility (|23|) on the critical isochore tp — for p = 0.5, u = 0.1 and three 
temperatures: above (r = —0.30810), below (r — —0.30811) and very close to the critical temperature (r = —03081037). 
obtained by numerical integration of Eq. (|22p . The dotted line is the critical fixed point (x* — —0.2848...) obtained by 
standard RG analysis of Eq. (J26J) . The dashed line is the low temperature fixed point x < — p — = —0.9142 . . . (see Ref. 

EH). 




Notice that the rescaled inverse compressibility x (|2H)) is simply expressed in terms of T t by x — The evolution 

equation for the "renormalized" quantities now acquires a typical RG structure: 



d_ i a_ 



where U(x) is given by Eq. (|22p. Equation (|26p admits a non trivial fixed point: i.e. a solution r*(£) independent 
of t, regular and not identically zero on the whole real axis £ G (—00,00). A fixed point corresponds to a diverging 
compressibility at ip = 0, leading to x* — ^s-|c=o ~ —0.2848... (for p — 0.5), which perfectly agrees with the 
position of the plateau in the evolution shown in Fig. [2] for r ~ r c , and therefore identifies the critical point. 
Linearizing the evolution equation (|26[) around r*(£) we can investigate the stability properties of the fixed point: 
letting r t (C) = r*(£) + (5r t (£) and keeping terms up to first order in ST t we find solutions of the form 5T t (() — e At 7(C) 
where A is the eigenvalue and 7(C) is the corresponding cigcnfunction. As usual in this context we find two unstable 
eigenmodes, characterized by positive eigenvalues [fjj governing displacements from the critical point along the critical 
isotherm and isochore. The first is odd in £ and can be expressed as 7(C) = with A G = 1/2, while the other, even 
in C, must be computed numerically. The RG structure of the HRT-OZ evolution equation (|2"rj|) shows that scaling 
and hyperscaling are correctly reproduced by the theory. Due to the known relation between eigenvalues and critical 
exponents Q 

S = i + A 

A 

■< - I < 27 > 

we find 5 = 5, as implied by scaling when the correlation exponent r\ vanishes, while the other critical exponents 
follow from the numerical determination of A e and are shown in Table I for few choices of the parameter p: HRT-OZ 
displays non classical critical exponents weakly dependent on p. It is in fact apparent that, contrary to the sharp 
cut-off HRT-OZ approach, here the asymptotic equation (|22[) does not acquire a universal form but retains memory of 
the microscopic model through the value of the parameter p. Therefore the critical properties of the three dimensional 
fluid within our approximation will generally violate universality and both critical exponents and scaling laws will 
depend on the precise value of p. This situation closely resembles the case of non perturbative approximations to 
the RG flow equations, whose results do depend on the cut-off function [T8|. Note however that, as shown in Ref. 
fl2L [T3|. the critical exponents are correct to first order in a e = 4 — d expansion of the HRT-OZ equation (j4|) and, 
within a systematic e— expansion of the full HRT hierarchy, the dependence on p disappears and universality is indeed 
recovered. The critical exponents can be also evaluated by numerical integration of the asymptotic evolution equation 
(1221) close to the critical point with no reference to any RG procedure. The results, shown in Fig. [3] for the case 
p = 0.5 and u = 0.1 conform to the expectation of the fixed point analysis previously outlined both above and below 
the critical temperature. 

The renormalized equation (|26[) also allows to determine other relevant properties of the critical regime. In Table 
I we show the compressibility amplitude ratio defined as t/jj = C+ / C_ where 

X {rMr)) = C± |Arp (28) 
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Exponent 







1 
1 


5 


71 
'1 


U 2 = C+/C- 


"Exact" 


0.110 


0.327 


1.237 


4.789 


0.036 


4.76 


HRT (p = 0.25) 


0.007 


0.332 


1.329 


5 





4.2 


HRT (p = 0.50) 


-0.003 


0.334 


1.335 


5 





4.5 


HRT (p = 0.75) 


-0.010 


0.335 


1.340 


5 





4.6 



TABLE I: HRT-OZ estimates of the critical exponents and compressibility amplitude ratio in three dimensions for few values 
of the parameter p < 1 compared to the exact values [l{| obtained by extrapolation of high-temperature series expansions. 




Log 10 |Ar| Log 10 |Ar 



FIG. 3: Left panel: Log-Log plot of the inverse compressibility vs. reduced temperature for p = 0.5 and u = 0.1. Upper points 
are the HRT-OZ results below the critical temperature in the homogeneous phase, extrapolated to the coexistence curve. Lower 
points are on the critical isochore ip = above r c . Right panel: Log-Log plot of the coexistence curve (squares) and of the 
spinodal (triangles) for p = 0.5 and u — 0.1. Lines mark the asymptotic power law behavior with the critical exponents quoted 
in Table I: 7 = 1.335 and = 0.334. 



X is the isothermal compressibility and Ar = r — r c . The + sign corresponds to r > r c and tp(r) = 0, while the 
— sign refers to r < r c and if(r) along the coexistence curve. Note that in the numerical solution of the evolution 
equation, the value of the inverse compressibility on the binodal, i.e on the point of discontinuity, is not directly 
available and must be obtained by extrapolation of the results at neighboring mesh points. This limits the possibility 
to get arbitrarily close to the critical point and to obtain accurate estimates of the amplitude ratio C/ 2 . 
The scaling form of the equation of state of the model can be also obtained in the form: 

x ~ 1=Ml ~ lw (w^ (29) 

The numerical results for p — 0.5 are shown in Fig. [4] together with a polynomial approximation of the "exact" scaling 
function from Ref. fl9l |. Other choices of p do not modify appreciably the good agreement between HRT-OZ and RG 
results. 



B. Below the critical temperature 

As previously noticed (see also Ref. [l3<]), within the smooth cut-off HRT-OZ scheme, the inverse compressibility 
develops a singularity for t — > 00 at any r < r c : it attains a finite limit on the binodal while vanishes identically inside 
the two phase region. The way x _1 drops to zero below r c is governed by a low temperature fixed point: 



x< = lim e 2t 



a 2 * 



< t— >oo dtp 2 

defined by the divergence of U(x) in Eq. (|22[) and explicitly given by x< = p — 2^/2. The agreement between the 
formal analysis and the numerical results shown in Fig. [2] fully confirms this interpretation. 

It is interesting to examine more closely the HRT-OZ evolution of the inverse compressibility below the critical 
temperature in order to understand how the inclusion of long wavelength fluctuations drives such a behavior, which 
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FIG. 4: Scaling form of the equation of state W(x) as defined in Eq. (I29|l for p = 0.5. Symbols: HRT-OZ results obtained 
by numerical integration of Eq. (|22[) . Full line: polynomial approximant of the exact scaling function The metric factors 
fixing the units of r — r c and ip are determined so as to locate the coexistence curve at Ar = — |(/j| 1 ^' 3 and x _1 = M* _ 1 a t T a - 




FIG. 5: Snapshots of the inverse compressibility as a function of density at four values of the evolution parameter. From 
bottom to top: t = 0, t = 0.5, t — 1.75 and t — > oo. Panel a refers to r = —0.27 above the critical temperature r c , panel b to 
r — —0.34 below r c . 



perfectly matches the expectations for the case of a scalar order parameter. In Figs. EKEb we contrast the evolution 
of two typical isotherms above and below r c by showing the inverse compressibility for four values of the switching 
parameter t, as obtained by the numerical integration of Eq. (|22[) . At t = the initial condition is a smooth function 
of the density which reproduces the Landau- Ginzburg mean field structure. Below the mean field critical temperature 
(r = 0) it displays a large unstable region, corresponding to the presence of a Van der Waals loop. In general, the 
inclusion of density fluctuations, governed by the parameter t, leads to a gradual increase of the inverse compressibility 
thereby reducing the width of the unstable region. Above the true critical temperature this eventually leads to the 
complete suppression of the instability (see Fig. [5k). Instead, below r c , long wavelength fluctuations severely affect 
the shape of the isotherms: i) x i when negative, flattens as a function of density and approaches zero at large t; 
ii) The isotherm sharpens at the boundaries of the interval where < developing a discontinuity, which defines 
the binodal; in) At large length-scales the coexistence region widens again and in a small density interval close to the 
binodal the compressibility becomes positive before dropping to zero. The origin of the flattening and sharpening of 
the HRT-OZ isotherm has been analyzed in detail in Ref. [l3| where the smooth cut-off HRT-OZ has been studied 
for a pure <I> 4 field theory. As anticipated.the smooth cut-off prescription overcomes one of the main problems of 
the previous implementations of HRT-OZ [7j, providing finite compressibility at coexistence. In this framework, we 
may investigate whether it is possible to locate the boundary of the metastable phases inside the coexistence curve. 
The spinodal curve is not rigorously defined within equilibrium statistical mechanics, the free energy being a convex 
function in the whole phase diagram, but a special feature of the evolution shown in Fig. [5Jd suggests a way to 
discriminate between instability and metastability. In a metastable state we expect that density fluctuations on a 
large but finite range do not drive the system towards phase separation, i.e. x" 1 remains positive until t does not 
exceed a (large) crossover value t s . Only when density fluctuations on a macroscopic scale are included, does the 
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FIG. 6: Boundaries of the region defined by a negative compressibility as a function of the evolution parameter t for p — 0.5, 
u — 0.1. Curve a refers to r = —0.27 > r c , curve b corresponds to r = —0.34 < r c and curve c to r — —0.308105 < c . Full 
dots mark the coexistence region, empty dots the spinodal. The dashed lines are guides to the eye. The width of the unstable 
region has a minimum at the near-critical temperature at t = 8.5 and <p = 0.0024. 

first order transition take place. In contrast, in the truly unstable region x the inverse compressibility is negative 
throughout the evolution. The curves dividing the regions of positive and negative x~ l m the (ip,t) plane at different 
temperatures are shown in Fig. [6j together with the extent of the metastable region obtained via the criterion 
previously outlined. According to this definition, the shape of the spinodal line in a ip vs. Ar diagram appears to 
be characterized by the same critical exponent f} of the coexistence curve, as shown in Fig [3J The value t* of the 
evolution parameter t where the width of the unstable region has a minimum, diverges when the critical point is 
approached. This corresponds to the divergence of a characteristic length-scale R c -j- e , which may be identified as 
the critical droplet radius on the spinodal line, according to the standard droplet model picture of nucleation [2ll |. 
This divergence follows a power law consistent with the scaling R c £ 4- |Ar| _ly with v — 7/2 (see Table I). 



Here we show some results obtained by the numerical integration of the full smooth cut-off HRT-OZ equation as 
derived in Appendix A . We systematically investigated the HCYF at three values of the inverse range parameter 
z: z — 1.8, which mimics the usual Lennard- Jones potential, z — 4 and z — 7 which are more appropriate for modeling 
colloidal suspensions. The numerical solution of the partial differential equation has been performed by a full implicit 
predictor-corrector finite difference method on a density mesh of 2000 points. The estimated numerical error cannot 
be appreciated on the scale of the figures. 



The coexistence curve can be immediately read-off the numerical solution of the equation, for in the two phase 
region the inverse compressibility vanishes. The spinodal curve instead is evaluated following the criterion discussed 
in Section (IV-B) and illustrated in Fig. [6l The results are shown in Figs. I7I8I9I together with the coexistence curves 
obtained by Monte Carlo simulations and SCOZA. 

The HRT-OZ coexistence curve agrees very well with the Monte Carlo data for z = 1.8 and z = 4 while some 
deviation appears at z = 7 where the MSA closure, which underlies the present implementation of the smooth cut-off 
HRT-OZ, is less accurate. As expected, the coexistence curve becomes flatter by reducing the range of the attractive 
interaction, as testified by the different scale of the temperature axis in the three cases. Note that the SCOZA results 
are always close to the HRT-OZ binodal, with some deviation in the critical region where the HRT-OZ coexistence 
line is flatter. Again, the largest discrepancy between HRT-OZ and SCOZA occurs at high z, however, even in this 
case, the two critical temperatures differ by less than 1%. 



V. 



NUMERICAL RESULTS AND PHASE DIAGRAMS 



A. Spinodal and binodal 
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FIG. 7: Coexistence curve and spinodal of the HCYF for z = 1.8: HRT-OZ (full circles), SCOZA (line) and Monte Carlo [T? 
(error bars). The HRT-OZ spinodal is also shown (crosses). 




FIG. 8: Coexistence curve and spinodal of the HCYF for z — 4. Notation as in Fig. [7] The Monte Carlo simulations are taken 
from Ref. [U. 




FIG. 9: Coexistence curve and spinodal of the HCYF for z — 7. Notation as in Fig. [7] The Monte Carlo simulations are taken 
from Ref. 0- 
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B. The liquid-solid phase boundary 



The fluid-solid phase equilibrium of the HCYF has been widely studied, mainly because it allows one to describe 
the qualitative change in the shape of the phase diagram brought about by changing the attraction range. As is 
well known, as the range is narrowed, i.e., as z is increased, one goes from the situation where the phase diagram 
has both a stable fluid- fluid critical point and a fluid-fluid-solid triple point, to that where the critical point moves 
into the fluid-solid coexistence region, thereby making the fluid-fluid transition metastable with respect to freezing, 
and there is no triple point [22L [23l |. The former case is found in atomic and molecular fluids, while the latter is 
characteristic of many colloidal systems, including protein solutions (24|. The disappearance of a stable fluid- fluid 
transition is in fact a general feature of short-ranged attractive interactions, but the Yukawa potential was, together 
with the Asakura-Oosawa one [2^], among the first to be considered in this respect. 

HRT is not a theory of freezing, but of course it can be used to predict the equilibrium fluid-solid phase diagram, 
if the free e nergy of the solid is provided by some different procedure. We have resorted to the simplest and most 
adopted one [25L [H, H3, Hl| , based on Barker and Henderson perturbation theory [2^| . This approach was originally 
developed for the fluid phase, but it can be straightforwardly employed in the solid one as well. According to it, the 
Hclmholtz free energy of a solid of particles interacting via a hard-sphere plus tail potential is given by: 



where A s is minus the excess free energy per unit volume of the solid divided by fcgT, is the corresponding 
quantity for the reference hard-sphere solid, gf(r) is the radial distribution function of the hard-sphere solid averaged 
over the solid angle, and w(r) is the potential tail, i.e., the attractive Yukawa in the present case. Following previous 
investigations [25|, [28| , Eq. (|30|) has been truncated at the second-order term A2 , which has been estimated by the 
so-called "macroscopic-compressibility" approximation, also suggested by Barker and Henderson [2^ |: 



where xf 1S the isothermal compressibility of the hard-sphere solid divided by the ideal-gas value. Comparison with 
Monte Carlo simulations on depletion potentials (3(| has shown that, although the estimate of A2 given by Eq. (|3~lj) 
deviates quite substantially from its exact value, which would involve knowledge of the three- and four-particle 
distribution function of the reference solid, the sum of A\ and A 2 given by Eqs. (|30I31|) is nevertheless remarkably 
close to the simulation results for the full free energy A s , unless the interaction is very short-ranged. In fact, the 
truncated expansion (|30f has proved to be more accurate in the solid than in the fluid 28], of course provided the 
solid in study has the same lattice structure as the reference one. 

Equations (|3QI3ip were already used together with the sharp cut-off version of HRT-OZ to study the fluid-solid 
equilibrium of non-additive hard-sphere mixtures, including the Asakura-Oosawa case [HJ]. They were also applied 
together with SCOZA to both the HCYF considered here [12] and a two- Yukawa fluid with competing attractive 
and repulsive interactions [33|. As in these studies, the Helmholtz free energy Af has been obtained by integrating 
with respect to density the equation of state of the hard-sphere solid given by Hall [13], starting from a density 
po = 0.736 p cp , where p cp — V2a~ 3 is the density at close packing, and setting the integration constant by the result 
—Af'/po — 5.91889 determined via numerical simulation by Frenkel, also reported in [32| . The radial distribution 
function of the hard-sphere solid gf(r) has been described by the parametrizations developed by Kincaid and Weis [Hj] 
starting from Weis' Monte Carlo simulations (26| and by Choi, Ree, and Ree [36] on the basis of their own simulation 
results. These parametrizations share the same functional form, and differ only for the specific dependence of the 
parameters on the average density. In both cases, we used the data for the neighbor distance and coordination number 
of the face-centered cubic lattice reported in [33] • The integrals in Eqs. (|30I31|) were truncated at r c = 8a, and a 
long-range correction was added by setting gff(r) = 1 for r > r c . Such a correction is actually irrelevant for the case 
discussed below. 

The fluid-solid phase boundary has been determined by equating the pressure and chemical potential of the solid 
given by Eqs. (|30I31[) to those of the fluid as predicted by the smooth cut-off HRT-OZ. As expected, we found that, 
as the inverse range z is increased, the freezing line becomes wider, and for large enough z it lies above the fluid-fluid 
coexistence curve. Figure ITUl shows the phase diagram for the case of marginal stability, such that the freezing line 
is tangent to the fluid- fluid coexistence curve at the critical point. This separates the stable fluid-fluid transition 
regime from the metastable one. We have reported the results obtained using both Kincaid and Weis [H| and Choi 
and Ree [36| radial distribution function gf (r). The two sets of results are nearly superimposed, save for some small 
deviations along the melting line. In both cases we obtain z = 5.6 as the value of marginal stability. This result is 
very close to that given by SCOZA, z = 5.7, via the same recipe used here to get Carnahan-Starling thermodynamics 




(30) 




(31) 
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FIG. 10: Fluid-solid phase diagram of the HCYF for z — 5.6. The fluid and solid phases have been described by HRT 
and perturbation theory respectively (see Eqs. (|30I31[0 . Full circles: HRT fluid-fluid coexistence curve. The asterisk marks 
the position of the fluid-fluid critical point. Squares: freezing and melting lines obtained using the hard-sphere solid radial 
distribution function gf{r) by Kincaid and Weis [3E|. Crosses: same as above using the gf{r) by Choi, Ree, and Ree [3^ |. 
Lines are a guide to the eye. 
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FIG. 11: Pressure divided by ksT as a function of density for z = 1.8 at three temperatures: T = 1, T — 1.5 and T — 2 (from 
bottom to top). Line: HRT-OZ results. Full points: Molecular Dynamics data form Ref. [39l ]. The two open dots mark the 
boundaries of the HRT-OZ coexistence curve at the lowest temperature, characterized by a flat portion of the isotherm. 

in the fluid phase, see Sec. Ill, and is similar to the value z = 6 determined by Hagen and Frenkel [22j via Gibbs 
Ensemble Monte Carlo simulations. In the case of SCOZA, using the hard-sphere direct correlation function Cfj(r) 
given by Waisman [38j ] instead of fixing its range a priori gave z — 6.05 [32], in closer agreement with the simulation 
result. It is likely that a similar trend would be shown by the smooth cut-off HRT-OZ with the same Cn(r). 

C. Equation of state and thermodynamics 

The equation of state of the HCYF can be easily obtained by integrating the HRT-OZ differential equation. The 
pressure divided by fc^T is compared to Molecular Dynamics simulations in Figs. QT] and [12] for z = 1.8 and 
z = 4 respectively. The agreement is very good showing that HRT-OZ is able to provide a consistent picture of 
thermodynamics throughout the phase diagram of the fluid. 

Finally, in Fig. [T3]we plot the specific heat as a function of temperature at the critical density for the three cases 
we have investigated z = 1.8, z = 4 and z = 7. The expected A-line can be clearly seen in the HRT-OZ data. Note 
the different growth of Cy at low temperature among the three cases. 
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FIG. 12: Pressure divided by ksT as a function of density for z — 4 at four temperatures: T = 0.7, T = 1, T = 1.5 and T = 2 
(from bottom to top). Line: HRT-OZ results. Full points: Molecular Dynamics data form Ref. [391 ]. 
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FIG. 13: Constant volume specific heat as a function of temperature at the critical density. From left to right z — 1.8, z = 4 
and z = 7. The critical temperature is identified by a vertical dashed line. 



D. Correlation functions 



Although the Hierarchical Reference Theory has been especially devised for investigating the thermodynamical 
properties of fluids, it also allows for the study of correlation functions through the link between thermodynamics and 
correlations implied by the closure relation discussed in Section III (|61 10[) . However, due to the MSA-like structure of 
the adopted scheme, we cannot expect an accuracy of the quality found for the thermodynamics. In Figs. 11411151 we 
compare the HRT-OZ results with Monte Carlo simulations for different values of z. The HRT-OZ result considerably 
underestimates the contact value of g(r), especially at large z or at low density, as a consequence of the "linear" 
relationship between c(r) and the potential outside the core radius. 

As a side comment, we remark that the structure of correlations in the two phase region is known to be analytically 
related to the correlation functions of the two pure phases: liquid (L) and vapor (V). According to the theory of first 
order phase transitions [4l| , the density correlations at coexistence are linear combinations of the density correlations 
for the two pure phases, leading to the following analytical structure of the radial distribution function at coexistence: 



9(r) = ( ^ 
P 



P~ Pv 

PL - PV 



[9L{r) - 1] 



pv\ 

p ) 



PL- P 
PL - PV 



[gv{r) - 1] + 



PiPL + Pv) - PLPV 

P 2 



(32) 



This shows that in the two phase region, the radial distribution function g(r) tends at large distances to a non 
trivial limit, which differs from the expected uncorrelated result g(r) — > 1. Such a behavior is in fact consistent 
with the divergence of the isothermal compressibility in the two phase region, via the compressibility sum rule: 
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FIG. 14: Radial distribution function of the HCYF for z = 1.8 at T = 1.5 and p = 0.4. Full line: HRT-OZ results. Points: 
Monte Carlo simulation from [171 ]. 
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FIG. 15: Radial distribution function of the HCYF for z = 4 at T = 0.5 and p = 0.816. Full line: HRT-OZ results. Points: 
Monte Carlo simulation from [ipl ]. 

ksTpXred = 1 + P J dr (g(r) — 1). Remarkably, the smooth cut-off HRT-OZ is able to capture this peculiar structure 
of the two body correlation. In fact, it can be shown [l4[ that the asymptotic limit of g(r) predicted by HRT-OZ 
precisely reproduces the last term of Eq. (|32|) . However, further analysis shows that the closure adopted in this work 
is not able to describe the full structure of correlations in the two phase region expected on the basis of the exact 
result (|32p : a spurious long range decay in g(r) appears inside the coexistence curve. 

VI. CONCLUSIONS 

In this work we have analyzed a simple Ornstcin-Zcrnike closure to the Hierarchical Reference Theory of fluids in 
the smooth cut-off formulation for hard core Yukawa fluids. Among the main advantages of this approach we can 
recall a consistent representation of the critical point, as customary in the HRT scheme, including non classical critical 
exponents and scaling laws, together with a correct description of the first order phase boundary. Contrary to the usual 
sharp cut-off HRT-OZ, here the inverse compressibility shows the expected discontinuity on the binodal which also 
allows to define the critical exponent 7' describing the divergence of the isothermal compressibility when approaching 
the critical point from the low temperature side. The expected long distance limit of two body correlations in the two 
phase region is also correctly reproduced by HRT-OZ. An estimate of the spinodal line can be extracted by studying 
the evolution of the inverse compressibility when density fluctuations of larger and larger wavelength are included. 
The thermodynamics of the HCYF agree well with simulations both near and far from the critical point: coexistence 
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curves and equations of state can be predicted with high accuracy by HRT-OZ. 

Despite such a remarkable performance of HRT-OZ, there is still room for significant improvement. Nowadays 
the most challenging applications of liquid state theory concern the so called complex fluids, often characterized by 
extremely short range attractive interactions and/or more structured potentials with both attractive and repulsive 
tails Q. The application of HRT to this class of systems requires further work: the closure to the exact HRT 
hierarchy considered up to now is clearly not adequate to describe complex fluids, for it assumes a linear relation 
between the space dependence of the direct correlation function and the two body interaction. Moreover the rather 
artificial splitting between a short range "reference" part of the potential vn(r) and a long range "perturbation" w(r), 
although customary in liquid state theory, limits the applicability of the method to potentials where repulsive and 
attractive contribution can be easily identified. The HRT-OZ description of correlations, already for moderately short 
range interactions, shown in Figs. [T4T[T"5l is not fully satisfactory and should be improved. In Section III we pointed 
out a natural way to generalize the adopted closure, which paves the way to more elaborate, flexible and general 
approximation schemes able to combine the virtues of the HRT approach and the accuracy of more sophisticated 
liquid state theories. We also mention a straightforward generalization of the smooth cut-off HRT method to binary 
mixtures, which is potentially very useful in the framework of complex fluids. 

From a more fundamental point of view, few problems still remain open, and should be addressed in future investi- 
gations. The HRT-OZ critical exponents are always affected by a systematic error (which is of about 10% for 7) clearly 
related to the adopted Ornstein-Zernikc structure of correlations which implies 77 = 0. This feature does not allow to 
apply the HRT-OZ approach in two dimensions. In order to go beyond this approximation, it is necessary to include, 
at least in some partial form, the effects of fluctuations on the momentum dependence of two body correlations: an 
investigation of approximation schemes to the second equation of the HRT hierarchy (O should be considered in the 
future. Hopefully, such a generalization will also provide a fully consistent picture of correlations in the two phase 
region. 

Finally, even within the present approximation scheme we did not address the subt le p roblem of the analytic 
structure of the free energy near the coexistence curve. According to the droplet model [21(, later confirmed by an 
exact result for the Ising model [42], an essential singularity is present at coexistence: the n th derivative of the free 
energy with respect to the chemical potential should scale as it)?/ 2 in three dimensions. Some preliminary numerical 
evidence suggests that an essential singularity is indeed present also within HRT-OZ but determining the precise 
scaling requires some further analytical and numerical effort. 

We are grateful to Nigel Wilding, Dino Costa and Giuseppe Pellicane for providing simulation results of the radial 
distribution function of the HCYF. 



Appendix A: The HRT-OZ equations and their numerical solution 



Here we discuss the formal structure of the smooth cut-off HRT-OZ equations (|4I6[) specialized to the three di- 
mensional (d = 3) Yukawa fluid (|11I12I13[) . We follow the notation of Ref. [43j (hereafter referred to as KJ) for the 
solution of the MSA equations. By substituting the form (JTTJ) into the evolution equation 0$ we get: 

^ = 2tt P 2 [Ci g t (s) + C 2 d s g t (s)] s=Z2 (33) 
Here gt(s) is the Laplace transform of r g t (r): 

g t (s) = / rg t (r)e- sr dr (34) 
Jo 

and the coefficients are defined as 

Ci = e -e- 2t (mt)-m) (35) 
C 2 = ±e z e- 3t m (36) 

where dots represent derivatives with respect to the evolution parameter t. According to Eq. (KJ-15) 



(37) 
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where 77 is the packing fraction 77 = np/Q. The functions t(s) and q(s) are given by Eqs. (KJ-16,KJ-17,KJ-18) in 
terms of (zx,z 2 ) and the parameters (a, b, (3x,f3 2 , dx,d 2 ) which satisfy the equations (KJ-10,KJ-11,KJ-12,KJ-13). This 
complicate set of six non-linear equations in six unknowns allows to find the MSA parameters as a function of the 
input variables: (p, z%, z 2 , K\, K 2 ). However, our specific problem is slightly different because K\ is actually unknown, 
depending on the consistency parameter At. Instead, we would like to express gt(s) in terms of {p, z, z 2 , K 2 ) and the 
compressibility which appears into Eq. ([6]). Luckily, the inverse compressibility is easily expressed in terms of the 
parameter a by 

p A' t = p J drc t (r) =1- a 2 (38) 

where a prime represents differentiation with respect to p. Therefore, in order to integrate the HRT-OZ differential 
equation (f3"3"]l we should solve, for each t, the five equations (KJ-10,KJ-11,KJ-12) and the second of the pair of equations 
(KJ-13) obtaining the five parameters (6, fix, f3 2 , dx, d 2 ) in terms of (p, zi, z 2 , A t , K 2 ). Then, by substituting these 
expressions into Eq. (|37)) we can evaluate the right hand side of Eq. (|33[) which becomes a partial differential equation 
(PDE) for the free energy density A t {p). Unfortunately, an analytic expression for the parameters (b, [3x , f3 2 , d\ , d 2 ) is 
not available and then we decided to adopt a different procedure. 

As a first step it is convenient to change the parametrization of the MSA equations because, during the HRT 
evolution, the inverse range of the potential z 2 — ze~ l becomes extremely small giving rise to singularities in the 
parameters (a, 6, j3t, (3 2 , dx, d 2 ). We first introduce /3, and defined by: 

k = M~ 2 (39) 
d,z\ = -1277 + A, (40) 

together with 

k 2 = K 2 z 2 A (41) 
which attains a finite limit for t — > 00. Then we note that the equations depend on the particular combination 

7 = 6- 1277^/3, 

i 

which will be used in place of b. The new set of parameters is then given by (a, 7, /3i, (3 2 , Ai, A2). Actually only the 
renormalization of f3 2 and A 2 are necessary but, for formal convenience, we adopted the same definition also for /3i 
and Ai. The Laplace transform of the radial distribution function is written as 

1 ( N(s) 



K.) = 3 + ^ (42) 



where the ideal gas term has been written explicitly. The key quantities N(s) and D(s) are expressed in terms of the 
parameters by the following set of definitions: 

S = E (- 12 ^ + A,) + fx(zi)} h-^- (43) 

^ Zi + S 

% 



r 



E^-^- (45) 

^- — ^ -4- _Q 



z, + s 

..2 



T (s) = a+{a + "f)s + s z T-12risY.o (46) 

q(s) = f 5 (s)r(s) (47) 

P = E-^-S r + ^ + *W (48) 

N(s) = r(s)/ 2 (-s)+r-a-7-sr + 1277E + 12v7P (49) 

D(s) = a + 7 s- I277 (s 2 P + sE ) (50) 

where we introduced the notation f n {x) for: 

n-X 



fn{x) = X- n 



3=0 J ' 
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which attains a finite limit as x — > 0. The two equations KJ-12 and the single equation KJ-13 now become, respectively: 



3 * - 

2 D(z 2 ) 



Z% - A 2 = 12-qzl 



k 2 z 2 = f3 2 D(z 2 ) 



while the first two equations KJ-10,KJ-11 are: 



1 - a 



12-q 



12r] 



f - 1 - E & A iM z i) + 12 ^E # *< M*) + E ^ 



% e 



(51) 

(52) 
(53) 

(54) 
(55) 



By writing explicitly the ideal gas contribution to the radial distribution function (|42p the evolution equation 
acquires an additional term which precisely generates the mean field part of the free energy. This term can be absorbed 
into the definition of a modified free energy density At : 



At = A t +A ld 



= A t + Ai, 



El 

2 

2ire 



Tz 2 



(56) 



where <p(k) denotes the Fourier transform of cj)(r) 
energy density is 



dA t 
dt 



2tt P 2 



Ci 



D(s) 



G 



(3w(r). Therefore, the evolution equation of the modified free 
D(s)d s N(s)-N(s)d s D(sY 



D( s y 



Now we define the new variable 



U t = 27T/9 2 



Tz 4 



fihM ti+W , lM 3 D(s)d s N(s)-N(s)d s D( S ) 
(2^{t) - ^(t)) z 2 — + z 2 



which basically coincides with the right hand side of the evolution equation (|33|) . which in fact simply becomes 



dA t 
dt 



27re 2 



Analogously, the compressibility equation ()38|1 is expressed in terms of the modified quantities as: 

- P A t = a 2 - pj^ z 2 e 2t T/j(t) 



(57) 



(58) 



(59) 



Combining Eqs. 

„2. 



and (|59p we finally get the evolution equation for the bare dimensionless inverse compressibility 



2aa — —pz 2 u'l 



(60) 



In summary, we start from the definition (J57J). We differentiate with respect to t, which is also contained in z 2 and in 
the six parameters (7, pi, A, ; , a). The derivative of the first five parameters are calculated by differentiating the five 
equations (|5H55p and solving the linear problem in terms of d. Finally, the derivative a is expressed in terms of u" 
via Eq. (|60p . This procedure provides a closed equation for ut which must be coupled to the five algebraic non-linear 
equations (|51ti55[) and to the definition (|57p in order to formally express the six parameters in terms of u t - An analytic 
solution of the set of equations (|5H55|57f is clearly impossible, but we can implement a Raphson-Newton iterative 
scheme at each step in t, which in fact converges very quickly. The resulting evolution equation for the independent 
variable Ut {p) has the general structure 



du t , . . . d 2 Ut 
= gi(u t ,t) + g 2 {ut,t)- s - T 



dt 



dp 2 



(61) 
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which is a quasi-linear PDE in two dimensions and it is stable provided (72 > 0. This condition must be checked 
during the numerical integration and depends on the specific definition of the cut-off function tp(t). In this work we 
have considered two different choices for the function ip(t): In the first part of the integration (0 < t < t*) we used 

4>(i) = (i + t/t y 2 (62) 

where t is a parameter. In the second part (t > t*) we have chosen 

ip(t) = (coshtr 2 (63) 

The transition point t* is obviously defined by the equation 1 + t*/to = cosht*. 

The initial condition at t = corresponds to a vanishing attractive tail Wt (r) and can be obtained by substituting 
the analytic solution of the MSA equation for a hard sphere liquid, equivalent to the Percus-Yevick (PY) solution, into 
the definition of u t (|57p . This form of the reference correlations is however rather inaccurate at high density, where 
the PY solution shows appreciable deviations from the Verlet-Weis form (44j. A better choice is to parametrize the 
hard sphere direct correlation function of the hard sphere liquid outside the core as a Yukawa tail. For convenience we 
decided to fix the inverse range of this additional Yukawa equal to z, i.e. the inverse range of the physical attractive 
part of the potential which defines the original problem. The amplitude of this additional contribution to the direct 
correlation function is fixed by requiring that the spatial integral of c(r) is consistent with the Carnahan-Starling 
equation of state [l5[ via the compressibility sum rule ([38]). This choice is in fact equivalent to setting the initial value 
of the parameter At to a density dependent non vanishing limit. The numerical solution of the PDE also requires two 
boundary conditions at the minimum and maximum density: at p = u t vanishes due to the definition (|57|l while at 
P = Pmax ~ 1 the MSA approximation, corresponding to A f = 0, provides a reasonable approximation. The PDE has 
been solved by a finite difference fully implicit predictor-corrector algorithm [45| with stability and accuracy checks 
during the numerical integration. The parameter t in the cut-off function has been chosen of order 1/z in order to 
preserve stability at all steps. 



Appendix B: The asymptotic form of the HRT-OZ equations for a Yukawa potential 



Let us first consider the simpler RPA based closure (fTT|) of the HRT-OZ evolution equation. By substituting Eq. 
(|17p into (|4]) and by rescaling the momentum variable q = fee*, we see that for t — > 00 the direct correlation function 
£t(fc) is evaluated in the k — > limit and can be therefore expanded: 

c t (k)^^~ bk 2 -^e^Uke 1 ) (64) 

where b attains a finite limit in the t — > 00 also at the critical point and is usually close to the range squared of the 
Yukawa tail while ipoo is given by Eq. (jSJ). Notice that even in the k — » limit, the direct correlation function does 
depend on the precise form of the two body interaction cj>(q) which, in our case is just the three dimensional Fourier 
transform of the Yukawa potential divided by kgT: cj>(q) = 4ir(3/(z 2 +q 2 ). After some algebra, the evolution equation 
for the modified free energy (|I8[) acquires the following form in three dimensions (d = 3): 



OAt e- 3t f d 3 q 24>(q)-4>'(q) 



Let us introduce the parameter 



(65) 



p = fez 4 /(4^oo) (66) 

measuring the ratio between the curvature of the potential and of the direct correlation function (times ipoo)- Then 
we rescale the inverse compressibility via Eq. (|20l) leading to 



= e -« *- r dy y2{1 + 2y2) (67) 



The integral is well defined provided the denominator does not vanish in the integration domain. Taking into account 
that p > we must have x + 1 > and either p + i>0orp + i<0 and x > p — 2^/p. The latter condition may be 
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satisfied only for p < 1. This integral is performed analytically by contour integration with the result (fl"9|) , valid for 
all p and x satisfying the above conditions. 

It is now natural to rescale the physical quantities (free energy and density) in order to get rid of the prefactors in 
Eq. (Ull): 



P- Pc 



-p- 1 / 2 ^t{ip) + 

1/2 



1 



-3i 



In terms of the new variables ^> t and ip we simply get 



(68) 



at 



d 2 ^ t 

dip 2 

" 3t U{x) 



where 



U(x) = 



^p + X + 2^/py / l + X 



(69) 



(70) 



which coincides with Eq. (|22|) . Differentiating the latter definition with respect to t, implicitly contained in x(t,ip), 
we have 



dU 
~dt 



dU 
dx 



2x 



dip 2 



(71) 



The equation is now written in quasi-linear form. We just have to formally obtain the coefficients by inverting the 
definition (|70|) and evaluating dU /dx. 



By defining the auxiliary quantity y = y/p(l + x) all the ingredients appearing in (I71| are given by: 



y --l 





P 












y = 




U 2 {\- 


P)- 




~P? 


u P -v 2 ^/u 2 - 


-P- 


-4-f 


-U 2 - 


dU 
dx 






D 


-3 






D = 




4- 


-P 









2U + ^Jp{U 2 - p + 4) 



(72) 
(73) 
(74) 
(75) 



The asymptotic analysis of the evolution equation with the full MSA-based closure (fl"U|) , which includes the core 
condition, is more involved. By inspection of the equations for the MSA parameters (|51ti55|) . it is natural to guess 
that 7, j3i, fa, Ai, A2 are all finite in that limit Actually, outside the coexistence curve, A2 — > 0, as shown by the 
numerical solution of the equations. In particular, in the critical region, we may assume that A2 — > and analyze the 
equations in the limit Z2 — > 0, while a/z2 and A2/Z2 are finite. To leading order in a and z%, the equations KJ-10, 
KJ-11 and the first of the KJ-12 just depend on (3i, Ai and 70 = 7 + ^2rj(3\. Therefore, the two equations (|52I53[) 
provide the limiting values of /?2 and A2/22: 



P 2 = (12//)- 1 -{a/z2 + 7o) + V(a/z2 + 7o) 2 + 24ryfc 2 



A 2 = 



10 ftAi/zi - 7o(l + ri/2) + /3i(12t7 - A 1 )/ 4 (z 1 ) + MM 
-12r) z 2 

^/{a/z 2 + 70) 2 + 2477/C2 



(76) 
(77) 



where &2 has been defined in (|4Tj) . When these expressions are substituted into the definition of ut (|57|) the singular 
contribution to u t becomes 



sinq 
U t OC Z2 



P2J0 - 4fc 2 



y/(a/z 2 +7o) 2 + 2477fc 2 
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besides an additive constant. Precisely at the critical point we expect that a ~ Z2, while inside the coexistence curve 
the ratio a/z 2 will tend to the value defined by the vanishing of the denominator: (a/z 2 + 70) 2 + 24rjk 2 = (recall 
that k 2 < 0) leading to a divergence of u s t ln9 / z 2 as t — > 00 and to a finite limiting value of A2. By comparing this 
expression with that obtained for a pure RPA (|T0|) , it is possible to identify the effective RPA parameters p and x as: 



To 2 



2477|fe 2 | 
a 2 

24t7|/c 2 |z 2 2 



(78) 
(79) 



The requirement < p < 1 then gives < 7g < 24ry|fc2| which should be verified in the region where phase transitions 
take place. 
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